load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "~/work/mylib/callvar.ncl" 

begin

  tp = "all" 
  tp = "tropic" 
  tp = "polar" 
  tp = "midlat" 
  tp = "CIRRUS2006" 

  fna = "K2009_"+tp+"_Ni.dat" 
  fnb = "K2009_"+tp+"_T.dat" 

  ni = asciiread(fna,-1,"float") 
  st = asciiread(fnb,-1,"float")

  ;;................................. 
  ;; unit conversion q
  ;;................................. 

  ni = ni * 1000. 

  nt = dimsizes(ni) 

  ni!0 = "time" 
  ni&time = ispan(1,nt,1) 

  st!0 = "time" 
  st&time = ispan(1,nt,1) 

  
  ni_210 = where(st.gt.205. .and. st.lt.215, ni, -999) 
  ni_210@_FillValue = -999

  ni_220 = where(st.gt.215. .and. st.lt.225, ni, -999) 
  ni_220@_FillValue = -999

  ni_230 = where(st.gt.225. .and. st.lt.235, ni, -999) 
  ni_230@_FillValue = -999

  ni_240 = where(st.gt.235. .and. st.lt.245, ni, -999) 
  ni_240@_FillValue = -999

  n210 = vnmiss_1D(ni_210) 
  n220 = vnmiss_1D(ni_220) 
  n230 = vnmiss_1D(ni_230) 
  n240 = vnmiss_1D(ni_240) 

;;...........................................................................
;; 210K 
;;...........................................................................
  qsort(n210)

  dm = dimsizes(n210)

  x10 = round(.10*dm,3)-1
  x25 = round(.25*dm,3)-1
  x75 = round(.75*dm,3)-1     ; at 0            
  x90 = round(.90*dm,3)-1

  pv = n210 
  pp = (/pv(x10),pv(x25),dim_median(pv),pv(x75),pv(x90),dim_avg(pv)/)

  print(" ")
  print(" ")
  print(" 205K percentile : ")
  print(" 10% : "+pp(0))
  print(" 25% : "+pp(1))
  print(" 50% : "+pp(2))
  print(" 75% : "+pp(3))
  print(" 90% : "+pp(4))
  delete(pv)

  pa = pp 

;;...........................................................................
;; 220K 
;;...........................................................................
  qsort(n220)

  dm = dimsizes(n220)

  x10 = round(.10*dm,3)-1
  x25 = round(.25*dm,3)-1
  x75 = round(.75*dm,3)-1     ; at 0
  x90 = round(.90*dm,3)-1

  pv = n220
  pp = (/pv(x10),pv(x25),dim_median(pv),pv(x75),pv(x90),dim_avg(pv)/)

  print(" ")
  print(" ")
  print(" 215K percentile : ")
  print(" 10% : "+pp(0))
  print(" 25% : "+pp(1))
  print(" 50% : "+pp(2))
  print(" 75% : "+pp(3))
  print(" 90% : "+pp(4))
  delete(pv)

  pb = pp 

;;...........................................................................
;; 230K 
;;...........................................................................
  qsort(n230)

  dm = dimsizes(n230)

  x10 = round(.10*dm,3)-1
  x25 = round(.25*dm,3)-1
  x75 = round(.75*dm,3)-1     ; at 0
  x90 = round(.90*dm,3)-1

  pv = n230
  pp = (/pv(x10),pv(x25),dim_median(pv),pv(x75),pv(x90),dim_avg(pv)/)

  print(" ")
  print(" ")
  print(" 225K percentile : ")
  print(" 10% : "+pp(0))
  print(" 25% : "+pp(1))
  print(" 50% : "+pp(2))
  print(" 75% : "+pp(3))
  print(" 90% : "+pp(4))
  delete(pv)

  pc = pp 

;;...........................................................................
;; 240K 
;;...........................................................................
  qsort(n240)

  dm = dimsizes(n240)

  x10 = round(.10*dm,3)-1
  x25 = round(.25*dm,3)-1
  x75 = round(.75*dm,3)-1     ; at 0
  x90 = round(.90*dm,3)-1

  pv = n240
  pp = (/pv(x10),pv(x25),dim_median(pv),pv(x75),pv(x90),dim_avg(pv)/)

  print(" ")
  print(" ")
  print(" 235K percentile : ")
  print(" 10% : "+pp(0))
  print(" 25% : "+pp(1))
  print(" 50% : "+pp(2))
  print(" 75% : "+pp(3))
  print(" 90% : "+pp(4))
  delete(pv)

  pd = pp 


;;...........................................................................
;; output data 
;;...........................................................................
  fno = "boxdata_hohn_obs.nc"

  system(" rm "+fno)

  flo = addfile(fno,"c")

  pa!0 = "pt"
  pb!0 = "pt"
  pc!0 = "pt"
  pd!0 = "pt"

  flo->BOX_210K=pa
  flo->BOX_220K=pb
  flo->BOX_230K=pc
  flo->BOX_240K=pd



end



